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Abstract A time-dependent model for the energy of a flaring solar active region 
is presented based on an existing stochastic jump-transition model ( jWheatland and Glukhov, 1998| | 
|Wheatland,"2008| fWEeatland, 2009] ) . The magnetic free energy of an active re- 
gion is assumed to vary in time due to a prescribed (deterministic) rate of energy 
input and prescribed (random) jumps downwards in energy due to flares. The ex- 
isting model reproduces observed flare statistics, in particular flare frequency-size 
and waiting-time distributions, but modeling presented to date has considered 
only the time-independent choices of constant energy input and constant flare 
transition rates with a power-law distribution in energy. These choices may 
be appropriate for a solar active region producing a constant mean rate of 
flares. However, many solar active regions exhibit time variation in their flare 
productivity, as exemplified by NOAA active region AR 11029, observed during 



October-November 2009 (Wheatland, 20101. Time variation is incorporated into 



the jump-transition model for two cases: 1. a step change in the rates of flare 
transitions; and 2. a step change in the rate of energy supply to the system. 
Analytic arguments are presented describing the qualitative behavior of the 
system in the two cases. In each case the system adjusts by shifting to a new 
stationary state over a relaxation time which is estimated analytically. The model 
exhibits flare-like event statistics. In each case the frequency-energy distribution 
is a power law for flare energies less than a time-dependent rollover set by the 
largest energy the system is likely to attain at a given time. The rollover is not 
observed if the mean free energy of the system is sufficiently large. For Case 1, the 
model exhibits a double exponential waiting-time distribution, corresponding to 
flaring at a constant mean rate during two intervals (before and after the step 
change) , if the average energy of the system is large. For Case 2 the waiting-time 
distribution is a simple exponential, again provided the average energy of the 
system is large. Monte Carlo simulations of Case 1 are presented which confirm 
the estimate for the relaxation time, and confirm the expected forms of the 
frequency-energy and waiting-time distributions. The simulation results provide 
a qualitative model for observed flare statistics in active region AR 11029. 
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1. Introduction 

Solar flares are magnetic explosions in the Sun's outer atmosphere, the corona, 
which occur intermittently. Local space weather effects due to large flares, in- 
cluding damage to electronics on expensive communications satellites, moti- 



vate attempts to predict flare occurrence (Odenwald, Green and Taylor, 2006 



|Committee On The Societal fc Economic Im pacts Of Severe S pace Weather Events, 2008[ ).| 
However, existing prediction methods are probabilistic, and not very reliable 
( | Wheatland, 2005} | Barnes and Leka, 2008] ). 

The flare mechanism is believed to be magnetic reconnection, which occurs 
at certain topological sites in the coronal magnetic field configuration of active 
regions around sunspots ( |Priest and Forbes, 2002) ). However, the reconnection 
process is incompletely understood, and detailed observations of individual flares 
reveal a diversity of effects which are often secondary phenomena, providing 
only indirect information on the underlying flare mechanism QBenz, 2008[ ). The 
statistics of solar flare occurrence have been examined in an attempt to provide 
insight into the flare phenomenon, and also to improve flare prediction. 

1.1. Observed flare statistics 

Two statistical distributions of interest are the flare frequency-size distribu- 
tion, and the waiting-time distribution. The frequency-size distribution Af(S) 
is the number of events per unit time and per unit size S, where 'size' refers 
to a measure of the magnitude of an event, for example the peak flux in X- 
ray wavelengths, or an estimate of the energy. The frequency-size distribu- 
tion is observed to be a featureless power-law over many orders of magnitude 
QAkabane, 1956| |Aschwanden, 2005[ ): 



M(s) = \ 1 ( 1 -i)sr 1 s-\ a) 

where 7 « 1.5-2 is the power-law index (the exact value depends on the choice 
of the measure for size) and Ai is the total mean rate of flaring for events 
larger than Si. Few exceptions to this simple power-law distribution have been 
reported [see however Wheatland! (|2010|> ] . although an upper rollover or cut-off 



is required on energetics grounds (Hudson, 1991 ). The waiting-time distribution 



is the distribution P(At) of times At between flare events. For individual active 
regions the waiting-time distribution is often observed to be consistent with a 
simple exponential: 

P(At) = X lC - XlAt . (2) 

where Ai is the mean rate of events of events defined by Equation ([T]). This model 
corresponds to flares occurring as independent random events at a constant mean 
rate, i.e. as a Poisson process in time ( |Moon et ai, 200"l) | Wheatland, 2001[ ). 
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However, the mean rate of flaring in active regions is often observed to vary, in 
which case a piecewise-constant Poisson model may be appropriate. The model 
waiting-time distribution is then a sum of exponentials corresponding to distinct 
intervals with different rates ([Wheatlan d and Litvinenko, 2002] ) : 

i 

where n±i = XuU is the number of events with size larger than S± corresponding 
to a rate Aij and an interval ti, and where Ni = nu is the total number of 
events. A power-law tail is observed in the combined waiting-time distribution 
for events from multiple active regions over longer periods of time on the Sun 
(Boffctta et ai, 19991, which also may be accounted for by the time-dependent 



Poisson model (Wheatland, 2000 Wheatland and Litvinenko, 20021, although 



some authors have argued for a non-Poisson interpretation ( Lepreti, Carbone, and Veltri, 200~lj ).[ 
Recently [Asclrwan den and McTiernanl (|2010j) demonstrated that a variety of 
observations for the longer-term waiting-time distribution are consistent with 
Poisson occurrence in time according to a simple functional form for the distri- 
bution of rates in time, which is a variant of the rate distribution presented in 
IWheatlandl(l2H0T))l . 

1.2. Flare statistics in two active regions 

To illustrate these ideas we consider two active regions: US National Oceanic and 
Atmospheric Administration (NOAA) active region AR 10486, from October- 
November 2003, and NOAA AR 11029, from October-November 2009. The data 
used are the soft X-ray event lists compiled by the US Space Weather Prediction 
Center (SWPCfJ The events are selected from whole-Sun 1-8 A flux measure- 
ments by the Geostationary Observational Environmental satellites (GOES). 
The peak flux of GOES events is routinely used to classify flares, and is the 
measure of size used here. The SWPC/GOES data are often used in studies of 
flare statistics, although they are less than ideal for this purpose because the lists 
are incomplete, and the peak fluxes in the lists are not background-subtracted 
[for a discussion see e.g. lWheatlandl(l200ip ]. However, the events for active region 
AR 11029 shown here are individually background subtracted, based on work in 



an earlier study Q Wheatland, 2010| ). 

Figure Q] illustrates the frequency- energy, and waiting-time distributions, for 
soft X-ray flares observed in active region AR 10486. This large and highly 
complex sunspot region appeared in the late stages of the maximum of the last 
solar cycle and produced a remarkable sequence of extremely large solar flares 
( |Dun et ai, 2007[ |Chumak, Zhang, and Guo, 2008[ ). Figure [1] shows the events 
in this region larger than peak flux Si = 3 x 10 -6 W m -2 (this choice should 
ensure the list of events is relatively complete, for events larger than Si). The 
region produced N\ = 43 events with peak flux larger than S±, including the 
biggest flare of the modern era (with listed peak flux 2.8 x 10~ 3 Wm~ 2 ), on 



1 Scc http: / /www. swpc. noaa.gov/ 1 
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4 November 2003. The upper panel shows the peak fluxes versus the recorded 
peak times for the events as a sequence of vertical lines, with a logarithmic 
scaling for flux, the 4 November flare being the tallest line. The middle panel 
plots the the number of events above size S versus size, which corresponds to the 
cumulative peak-flux distribution C(S), defined in terms of the frequency-peak 
flux distribution by 

C(S) = — / Af(S')dS', (4) 
i J s 

where T is the duration of the observing interval, and S > Si. The vertical line 
in this panel indicates Si- The events are approximately power-law distributed 
in peak flux (a straight line in this representation), consistent with Equation (JlJ. 
The lower panel plots the waiting-time distribution for the same events, again 
as a cumulative distribution i.e. the number of waiting times larger than a given 
time versus waiting time. This corresponds to the integral of -P(Ai) with respect 
to At. The distribution is presented with a log-linear scaling, and reveals that 
the events are approximately exponentially distributed (a straight line, in this 
representation), consistent with Equation ©. 

Figure [T] shows that AR 10486, although remarkable in terms of the size of 
particular events, was unremarkable statistically. It produced flares according to 
a power-law frequency-peak flux distribution at a constant mean rate during its 
transit of the disk. This suggests that the physical conditions underlying flaring 
are approximately time- independent. 

Figure [5] shows GOES event data for active region AR 11029, which ex- 
hibits time variation in its mean rate of flaring. This spatially small but highly 
flare-productive active region emerged on the disk in October 2009, during an 
extended interval of low solar activity. The statistics of events in this region were 
investigated in [Wheatland] (|2010|) . The presentation of Figure [5] is the same as 
FigurcQ] but in this case the data are individually background-subtracted, which 
is particularly important because the events are small. There are Ni = 56 events 
above a backgrounded-subtracted peak flux of Si = 10~ 7 Wm -2 , as shown in 
the upper panel. This panel suggests that the flaring rate is high for two days (26 
October and 27 October), and relatively low at other times, an interpretation 
supported by a Bayesian rate analysis ( |Scargle, 1998 Wheatland, 2010[ ). The 



changes in the rate are quite dramatic: a relatively sudden increase by a factor 
of about ten, and a sudden decrease by about the same factor. The lower panel 
plots the cumulative waiting-time distribution, which shows a double exponential 
form, consistent with Poisson occurrence at two different mean rates (a low and 
a high rate), as represented by Equation ([3]) with two intervals/rates. The middle 
panel in Figure [5] plots the cumulative peak-flux distribution for the flares in AR 
11029, and is suggestive of a rollover around 10~ 6 Wm -2 . IWheatlahdl f|2010[) 
applied Bayesian model comparison to show that a power-law plus upper rollover 
model is much more probable for this data set than a simple power-law model, 
and argued that the departure from a power law may reflect the finite storage 
of magnetic energy in this (small) region. 

Time variation in mean flaring rate is commonly observed in active regions 
d Wheatland, 2001] ), and is often associated, as in the case of active region AR 
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Figure 1. Flares in the GOES events lists for NOAA active region AR 10486 above peak soft 
X-ray flux 3 X 10 — 8 Wm" 2 (43 events). Upper panel: Schematic showing the peak flux of each 
event at the event time. Middle panel: The frequency-peak flux distribution, as a cumulative 
plot in a log-log representation. Lower panel: The waiting-time distribution as a cumulative 
plot in a log-linear representation. 



11029, with a change in the photospheric magnetic complexity of a region ([Wheatla nd, 2010] ) .[ 
On October 26 the region increased in photospheric magnetic complexity (be- 
coming a /3~7 region, in the Mt Wilson classification), coincident with the 
increase in flaring rate. Photospheric magnetic field changes are likely to be 
reflected in changes in the magnetic field configuration in the corona, which may 
facilitate or inhibit rcconncction, and hence change the flaring rate. It is plausible 
that the more interesting statistics observed for AR 11029, by comparison with 
AR 10486 (Figure [I) are related to the smaller physical smaller size of the 
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region, which may make the coronal magnetic field configuration more sensitive 
to change. 



Active region AR 1 1 029 
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Figure 2. Flares in the GOES events lists for NOAA active region AR 11029 with a back- 
ground-subtracted peak flux above 10 — 7 Wm -2 , a total of 56 events (Wheatland 2010). The 
presentation of the data is the same as in Figure 1. 



1.3. Models for flare statistics 

Models for solar flare statistics attempt to account for the observed distributions. 
Most of the models are motivated either by accounting for energy balance in an 
active region ( |Rosner and Vaiana, 1978| |Litvinenko, 1994| |Craig, 2001[ ) or by ex- 
plaining the power law in the frequency-size distribution using self-organized crit- 
icality/cellular automata ('avalanche' models) QLu and Hamilton, 19 91; Cha rbonneau et ai, 200l} [ 
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|Hughes et ai, 2003[ ), although the two pictures are not mutually exclusive. The 
avalanche model produces a power-law frequency-size distribution below an up- 
per rollover set by the size of the cellular automata ( Lu et al., 1993[ ) , and has an 



exponential (Poisson) waiting-time distribution, if the system is subject to a con- 
stant mean rate of driving ([Biesccker, 1994; Wheatland, Sturrock, and McTiernan, 1998] ) .[ 
but may produce other waiting-time distributions with time-dependent driv- 



ing ( Norman et ai, 2001 1. Energy balance models are designed to explain the 
power-law frequency- energy distribution, and generally assume Poisson flare oc- 
currence, but simple models in which active region energy is completely removed 
by a flare imply a relationship between waiting times and energy, which is not 
observed ( |Lu, 1995] ) . A generalized energy-balance formalism (a 'jump-transition 
model') which avoids this problem was introduced bv lWheatland and Glukhovl 
flTMEl), and further developed in IWheatlahdl ([2U05)) and IWheatlandl ([2009")) . The 
jump-transition model is applied in this paper. 



The general jump-transition model (Wheatland and Glukhov, 1998 Wlicatland, 2008 | 



|Wheatland, 2009[ ) describes the magnetic free energy E = E(t) of an active 
region at time t in terms of secular energy input at a rate /3(E, t), and stochastic 
transitions downwards in energy (flares) at a rate a(E,E' ,t) per unit time 
and per unit energy for transitions from E to E' . The model may be pre- 
sented either in terms of an intcgro-partial differential equation describing the 
probability distribution P(E,t) for the energy (Wheatla nd and Glukhov, 1998| 



Wheatland, 2008 ), called a 'master equation' (van Kampen, 1992 Gardiner, 2004),| 



or in terms of a stochastic differential equation describing the time evolution of 
E = E(t) ( |Wheatland, 20091 ). 

Modeling using the jump-transition formalism to date has considered only 
time-independent cases, namely a time- and energy-independent rate of energy 
input: 

P{E) = ft, (5) 

where B is a constant, and time- independent power- law distributed flare tran- 
sitions!^: 

a(E, E') = a (E - E'ye(E - E' - E c ), (6) 

where ao is a constant, E c is a low-energy cutoff, 9 (x) denotes the step function, 
and 7 = 1.5 is the observed power-law index in the flare frequency-energy dis- 
tribution ( |Aschwanden, 2005[ ). The model flare frequency-energy distribution is 
given by 

/>oo 

Af(E) = / P(E')a(E' , E' — E)dE' 
Je 

/>oo 

= a Q E-~< / P(E')dE' for E > E c , (7) 
Je 



z Wheatland (2008) and Wheatland ( 2009) considered also the case with an additional factor 
E s multiplying the transition rates given by Equation JHJ, but the results suggested that 
Equation (0 is the preferred model. 
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and is zero for E < E c . Equation ([7]) describes a power law with index 7 above 
E c and below an upper rollover set by the largest energy the system is likely 
to attain, which has an approximate lower bound given by the estimate for the 
average or mean energy of the system ( [Wheatland, 2008[ ): 

1/(2-7) 

(8) 



For a sufficiently large rate of energy supply with respect to the rate of flare 
transitions, we have 

a E c 

and hence £ 3> E c , i.e. the system mean energy is large. In that case most 
flares do not significantly deplete the free energy, and then, according to Equa- 
tion ([7]) , the observed flare frequency- energy distribution is a simple power law. 
The rollover is not relevant unless a very long time history of flaring (including 
infrequent, very large events) is observed. In the case £ 3> E c also, flares occur 
as a Poisson process in time, i.e. the waiting-time distribution is a simple expo- 
nential. This may be understood by noting that the total flaring rate assuming 
the system has energy E is 

X(E) = f a(E,E')dE' 
Jo 

= a (^ 7+1 -^ 7+1 )/(7-l) for E>E C , (10) 

and X(E) = for E < E c . When £ ^> E c we can assume E ^> E c , in which case 
Equation (fTU|) is well approximated by 

\(E) = -^-E~^ +x . (11) 
7-1 

The right-hand side of Equation (|llj) is energy- and hence time-independent, so 
flares occur as a simple Poisson process in time. However, if flares deplete the 
energy of the system, then the system energy may become comparable to E c , 
in which case the total mean flaring rate depends on the energy of the system 
and varies in time, and the waiting-time distribution departs from a simple 
exponential (Wheatland, 2008; | Wheatland, 2009j ) . These comments show how 
observed flare statistics provide insight into magnetic energy balance in an active 
region, in the context of a model. 

The jump-transition model is a general formalism describing time-dependent 
as well as stationary situations. In this paper we apply the model to time- 
dependent cases for the first time, to investigate how time variation affects 
the model event statistics, in particular the frequency- energy and waiting-time 
distributions. We focus on two cases: a sudden change in the flaring rate; and a 
sudden change in the energy supply rate. The first case may be appropriate to 
describe, e.g., active region AR 11029 on 26 October 2009 (see the upper panel 
of Figure ■ The second case may be appropriate to describe an active region 
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which starts to grow as a result of the emergence of new magnetic flux, but which 
retains a given magnetic configuration and does not change its flaring rate. 

The layout of the paper is as follows. Section [5] briefly re-iterates the details 
of the jump-transition formalism (Section I2.1[) . describes the flare-like choices 
for modeling time- variation considered here (Section 12. 2p . and then presents 
analytic arguments which allow the general behavior of the system to be deduced 
(Section l2.3[) . Section [3] presents Monte Carlo simulations of the time-dependent 
system for the case of a sudden increase in the flaring rate, which provide 
a qualitative model for the observed behavior of active region AR 11029. A 
brief description of the numerical methods is given ( Section 13. 1|) . followed by an 
account of the results (Section 13. 2[) . Section [4] presents conclusions. 



2. Model 

2.1. General master equation and stochastic differential equation models 

I Wheatland and Glukhovl (TIM5| . IWheatlandl ([2005)1 . and IWheatlandl ([MIS]) de- 
veloped a stochastic jump-transition model for the free magnetic energy of a 
solar active region, which is described by a master equation, or an equivalent 
stochastic differential equation ( |van Kampen, 19921 |Gardiner, 20040 . 

In the master equation approach, the probability distribution P(E,t) for the 
energy E of the system at time t is given by the solution of the integro-partial 
differential equation 

= --^l(3(E,t)P(E,t)}-\(E,t)P(E,t) 

+ / P(E',t)a(E',E,t)dE', (12) 
Je 

where j3(E, t) is the rate of energy input to the system, a(E, E', t) is the rate of 
flare jumps in energy from E to E' per unit energy, and 



X(E,t) = ( a(E,E',t)dE' 
Jo 

is the total rate of flaring when the system has energy E. If the energy supply 
rate and the transition rates do not vary with time then the time-independent 
version of Equation ([12")) applies (obtained by setting d/dt = 0). 

The observable distributions of interest are the frequency-energy distribution 
and the waiting-time distribution. The model frequency-energy distribution is 

/•OO 

Af(E,t)= P{E',t)a{E',E' - E,t)dE'. (14) 
Je 

In the time-independent case, the model waiting-time distribution may be ex- 
pressed in terms of the stationary probability distribution P{E) for the system 
energy and the solution to an auxiliary first-order partial differential equation, as 
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shown by |Daly and Porporato| (|2007p in the context of general jump-transition 
modeling. The details of this procedure are given in IWheatlandl ([2008) and 
IWheatlandl (|2009|) but are omitted here. In the general time-dependent case 
there does not appear to be a straightforward way to obtain the waiting-time 
distribution from the master equation approach. 

IWheatland and Glukhovl (|T??9"gf and IWheatlandl pTOS)) solved Equation dHJ) 
for the time- independent flare-like choices discussed in Section [TJ namely a 
constant rate of energy supply and power-law distributed transition rates [Equa- 
tions §5§ and ©] , and the frequency-energy and waiting-time distributions were 
investigated for these choices. An efficient numerical method for solving the 
time- independent master equation was given in [Wheatland (2008). 

Moments of the master equation, i.e. averages over energy, provide insight 



into the general behavior of the system (Wheatland and Litvinenko, 20011. The 



first moment, obtained by multiplying Equation (|12j) by E and integrating with 
respect to E, gives 

±(E) = {P)-(r), (15) 



(16) 



df 

where 

/>oo 

(E) = / EP(E,t)dE 
Jo 

is the mean energy, 

POO 

(l3) = / P(E,t)P(E,t)dE (17) 
Jo 

is the mean rate of energy supply, and 

/•OO 

(r) = / r(E,t)P(E,t)dE (18) 



(19) 



Jo 

is the mean total rate of loss of energy to flaring, where 



r(E, t)= ( (E- E') a(E, E', t)dE' 
Jo 



is the rate of loss of energy due to all jumps. Equation (fTS"]) describes how the 
mean energy of the system changes in response to time variation in the rate of 
energy supply or in the rates of flare transitions. 

The equivalent stochastic differential equation to Equation (TT2"j) is ( |Daly and Porporato, 200"7| )[ 



dt 
where 



dE 

P(E,t) - A(E,t), (20) 



N(t) 

A(E,t) = J2 AE i S ( t - t i) ( 21 ) 
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is the total loss in energy due to flaring up to time t, with S(t) being the delta 
function, N(t) being the number of events which have occurred up to time t, and 
where the event times ti are defined by the Poisson process with time-dependent 
rate A = A [E(t), t]. The factors AEi are the jumps downwards in energy at each 
flare, which follow the distribution h(AE, E,t), defined by 

a(E, E - AE, t) = X(E, t)h(AE, E, t), (22) 

and satisfying the normalization condition 

E 

h(AE,E,t)d(AE) = 1. (23) 

The stochastic differential equation (|2T))) is simulated for the time-independent 
flare-like choices, using an efficient Monte Carlo method, in [Wheatland (12009|) . 

2.2. Choices for time variation 



As a simple time-dependent generalization of the flare-like model defined by 
Equations ([5]) and (O, we consider the choices 

a(E, E', t) = a (t)(E - E'y i d{E - E' - E c ) (24) 

with a power-law index 7 = 1.5, and 

A, (t). (25) 

Hence we consider time-modulated flare transition rates with the same power- 
law functional form for the change in energy previously considered, and a strictly 
time-dependent energy supply rate. We further restrict attention to the following 
two cases. 

Case 1: a step change in the transition-rate coefficient ao(t), 

a (t) = aoi + (a 2-aoi)6(t-T), (26) 

with no change in the energy-supply rate coefficient Po(t) = /3qi. 
Case 2: a step change in the energy-supply rate coefficient, 

A>(t) = A>i + (A>a-A)i)0(*-T), (27) 

with no change in the flare transition-rate coefficient ao(t) = ao\. 

In Equations ((26)) and (|27p the factors «oi and Am (with i = 1,2) are constants, 
and t = T denotes the time of the step change. 

Figure [3] illustrates the two cases, with the upper row showing Case 1, and 
the lower row Case 2. Case 1 may correspond physically to a change in the 
coronal magnetic configuration, which permits enhanced reconnection rates to 
occur, and Case 2 may correspond to increased sub-photospheric driving, e.g. 
the emergence of new magnetic flux. Case 1 provides a simple model which may 
account for the observed behavior of active region AR 11029 on October 26 (see 
Section [TTU and Figure 0). 
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Figure 3. The two simple cases of time variation in the flare model considered. The upper 
row is Case 1, a step variation in the flare transition-rate coefficient cto(t), and no change in 
the energy supply rate /3o(*)i an d the lower row is Case 2, a step change in the energy-supply 
rate coefficient, and no change in the flare transition rate coefficient. The parameters ri and 
T2 are the relaxation times, described in Scction l2.3l 



2.3. Analytic considerations 

For both Case 1 and Case 2 we assume that the system is in a steady state 
prior to the time t = T of the step change. In other words we assume the 
system has had constant flare transition and energy supply rates (defined by the 
coefficients «oi and /?oi respectively) for an extended interval prior to t = T. 
The energy distribution Pi (E) describing the steady state is the solution to the 
time-independent master equation with these parameters. The average energy 
over this distribution is given approximately by Equation (J5J , namely 

Vaoi/poi/ 

In both cases the system is in a non-steady state immediately after the change. 
In the context of the master equation description, this means that the energy 
distribution is inconsistent with the solution to the time-independent master 
equation for the new (constant) rate coefficients. In the context of the stochastic 
differential equation description, it means that the energy is not a typical energy 
for the system. The evolution of the system in the non-steady state is described 
by a time-dependent energy distribution P(E, t) which is a solution to the time- 
dependent master equation (|12j) . or by a specific energy trajectory E = E(t) 
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obtained by solving the stochastic differential equation (j2TJ)) . A new steady 
state is eventually achieved in each case, characterized by new distributions 
P2i(E) (for Case 1) and P22(E) (Case 2), which are solutions to the time- 
independent master equation with parameters defined by rate coefficients ctQ2, 
/3oi, and aoi, A)2j respectively. The mean energies of these distributions are given 
approximately by 



£21 = 



( 2-7 



1/(2-7) 



and £22 = 



2-7 



1/(2-7) 



(29) 



respectively. The characteristic times for achieving a steady state are referred to 
as 'relaxation' times, and are labeled t\ and T2 for Cases 1 and 2 respectively. 
These times are indicated schematically in Figure [3] by the vertical dashed lines. 
The interval T < t < T + Ti (with i = 1, 2) is the relaxation interval in each case, 
during which time the mean or peak of the energy distribution P(E, t) shifts to 
a new value. For the specific examples shown in Figure O i.e. an increase in the 
flare transition rates for Case 1, and an increase in the energy-supply rate for 
Case 2, the peak of the energy distribution shifts to a lower energy, and to a 
higher energy, respectively. The locations of the peaks are defined approximately 
by Equations (|28p and (|2"9")l . Figure 2] illustrates these changes, showing the initial 
steady-state distribution Pi (E) and the final steady-state distributions P21 (E) , 
and P22(E), and the values £1, £21, and £22 - This is a schematic diagram, not the 
result of a calculation, but the distributions have been drawn to approximately 
match the functional forms observed for numerical solutions of the steady-state 
master equation ( |Wheatland and Glukhov, 1 998| |Wheatland, 2008[ ). 



P(E) 




£21 



£1 



£ 22 log E 



Figure 4. A schematic showing the steady state distribution Pi(E) before the step change 
in parameters, and the steady state distributions after relaxation: P2i(E) (for Case 1), and 
Pl2(E) (Case 2). These specific examples correspond to Figure [3] i.e. an increase in the flare 
transition rates for Case 1, and an increase in the energy supply rate for Case 2. 

The relaxation times for the two cases may be estimated based on the first 
moment of the master equation, presented in Section [2.f I Substituting the time- 
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dependent flare-like choices Equations and (f2~5j) into Equations 
gives the more specific form for the moment equation 



dt 



(E)=Po(t)-2ao(t)(E)" 



(30) 



where we assume E >• E c and (E 2 ~ J ) m {E) 2 ~~* . Equation ([3U|) is a nonlin- 
ear ordinary differential equation. Implicit analytic forms for solutions may be 
written down for the specific cases of interest (7 = 1.5 and time-independent 
rate coefficients), but it is simplest to estimate a relaxation time directly from 
Equation ([3U)l . 

For Case 1, the system shifts from the approximate average energy £\ to 
the approximate average energy £21 in the relaxation time t\ with parameters 
cto(t) = ao2 and /?o (t) — P01 during the interval of evolution, so Equation (130[) 
implies the approximate relationship 



|g 2 i - £1 

n 



An - 2a 2^i 



7-1 



Substituting Equations (|28p and (|29[) into Equation ([3~T]) leads to 



71 



A)( a oi + Q02) 



(31) 



(32) 



assuming ao2 7^ ctoi and using 7 = 1.5. 

Similarly, for Case 2, the system shifts from the approximate average energy 
£1 to the approximate average energy £22 in time T2 with parameters ao(t) = aoi 
and /?o(t) = P02, so Equation (|30p implies 



\£ 



22 



T"2 



which leads to 



T-2 



/3 02 - 2a i^i 



An + /3q2 
An 2 ; 



7-1 



(33) 



(34) 



assuming A)2 ^ An and using 7 = 1.5. Equations (|3"2")l and (jM)) are accurate if 
the changes are small. If the mean energy of the system increases substantially 
in each case the expressions give overestimates because of the use of the initial 
small mean energy to replace the time- varying average energy (E) on the right- 
hand side of Equation ([3TJ]) . Similarly the expressions give underestimates if the 
mean energy decreases substantially. 

We can also deduce general analytic results for the model waiting-time and 
frequency-energy distributions. Substituting the time-dependent flare-like choice 
Equation into the general expressions for the total flaring rate, Equa- 

tion (fT3"|) . and for the flare frequency- energy distribution, Equation (fl"4")l . gives 



X(E, t) = ao(t) {E^ +1 E^ +1 ) /( 7 - 1) 



(35) 
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and 

poo 

M(E,t) = a (t)E-~< / P(E',t)dE' (36) 

J E 

respectively. These expressions apply for E > E c : the total rate and the frequency-^ 
energy distribution are both zero if E < E c . 

The frequency-energy distribution defined by Equation (|3"6"|) is a power law 
with index 7 for flare energies less than a timc-dcpcndcnt rollover, defined by 
the largest energy the system is likely to attain at a given time. Prior to the 
step changes considered here, the rollover has the approximate lower bound S\. 
After the step changes the lower bound is £21 (for Case 1) and £22 (Case 2). 
If the active region has a very large mean energy before and after the change, 
then most flares involve decreases in energy substantially less than the mean 
energy, and the rollover is not observed in a short interval of observation. The 
frequency-energy distribution is then given approximately by a simple power 
law: 

JVoi(E) = a 0i E-i, (37) 

before (i — 1) and after (i = 2) the step change. 

The waiting-time distribution given by Equation (|35[) is in general the waiting- 
time distribution for a time-dependent Poisson process with a rate A = X[E(t), t\, 
which has a well-known form ( |Wheatland and Litvinenko, 2002 ). If the sys- 
tem has a very large mean energy, so that E 3> E Cl Equation (|35p may be 
approximated by 

X(E,t) =a (t)£ r : 7+ 7(7-l), (38) 

and time dependence enters only through the coefficient ao(t). In this case the 
total flaring rate has the approximate constant values 

A01 = a i-E- 7+ 7(7 - 1) and A 02 - a m E^ +1 / ^ - I) (39) 

before and after the step change respectively, so the waiting-time distribution is 
a simple exponential before and after the change, corresponding to Equation ([2]) 
(the waiting-time distribution for a time- independent Poisson process). Before 
the step change the exponent in the exponential is A01, and after it is Ao2- 
The waiting-time distribution constructed for events both before and after the 
step change is in general a double exponential, specified by Equation §3§ with 
two intervals and rates. For Case 2, the waiting-time distribution has the same 
exponential form (with exponent A01) before and after the step change, because 
the flare transition rate coefficient does not change. If the system mean energy 
is small, the distributions may depart from exponential forms, as discussed in 
Section [1.31 

3. Stochastic modeling 

To illustrate the time-dependent model, and to confirm the qualitative ana- 
lytic results given above, we present Monte Carlo solutions to the stochastic 
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differential equation formulation [Equations ([2T))) - ([2"3"|) ]. following the approach 
presented in [Wheatland (2009). We consider only Case 1, since it exhibits the 
more interesting event statistics, and provides a simple model for the observed 
behavior of active region AR 11029 (Section 11.21) . 

3.1. Numerical method 

Full details of the Monte Carlo method are given in [Wheatland (20M[). and here 
we present only a brief description appropriate for the step change modeling of 
Case 1. 

A simulation is started with an energy E s at time t s . The system is assumed 
to have a constant flare transition-rate coefficient aoi and a constant energy- 
supply rate coefficient /3oi- Prior to the hrst flare the system evolves in energy 
according to Equation (j2"0")) without the loss term, so the energy as a function of 
time is 

E*(t) =£?.+ An (*-*.)• (4°) 

The total (expected) flaring rate during this time is A* = X[E*(t),t], where 
X(E,t) is given by Equation (|35[) with otait) = aoi- A random Poisson waiting 



time At is generated corresponding to this rate ( Wheatland and Craig, 2006 
Wheatland, 2009]) , which defines the end time t e = t s + At when a jump transi- 
tion occurs. The energy prior to the jump is E e = E*(t e ). A random jump of size 
AE is generated from the distribution h(AE, E e ,t e ) defined by Equation (j22|) . 
as explained in [Wheatland] (|2009p . Once AE is calculated, the whole process 
is repeated, with the new starting time t s = t e and the new starting energy 
E s = E e — AE. The process is repeated n\ times to give a time history of energy 
over this number of jump transitions for the interval prior to the step change. 
The time of the last transition, which we label t\, defines the time T of the step 
change. The whole process is then repeated again ni times, with the new flare 
transition-rate coefficient ao2, to give a time history for the interval ti after the 
step change. The starting energy for the system following the step change is the 
energy after the last jump transition prior to the step change. 

This process represents a single simulation. Ensemble averages over repeated 
simulations allow comparison with the master equation formulation of the model, 
and the moment equation. 

For the purposes of the simulations we introduce non-dimensional parameters 

"TrT E — t — a 0t t s — f3 01 t s 

E =TT> < = "T' a ° l = ~> &>1 = —FT- ( 41 ) 

where t s is an arbitrary scale time, and the choice 7 = 1.5 is explicitly shown. 
We choose t s = E c //3qi without loss of generality, so that j3 01 = 1. The non- 
dimensional mean energies for Case 1 are then 

£1 = 1/oq! and £ 12 = l/aj^, (42) 

and the estimate for the relaxation time is 

r i ~ , 2 ■ ( 43 ) 

4a ia 02 
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3.2. Results 

A simulation of Case 1 is performed. The chosen parameters are cfoi = 5 x 10~ 3 
and ctQ2 = 5 x 10~ 2 [also /3 01 = 1, as explained following Equation (|41|]. These 
values are chosen to mimic solar active region AR 11029 ([Wheatla nd, 2010[ ), 
which exhibited an approximate ten-fold increase in its flaring rate on 26 Octo- 
ber 2009 (see Section [TT21 and Figure [2]). The corresponding approximate mean 
energies are £\ = 10 4 and £2 = 10 2 , the approximate total flaring rates [given 
by Equation (|39|) ] arc A01 = 25 i = 10~ 2 and A02 = 2a 2 = 10 -1 , and the 
approximate relaxation time is f\ = 1.1 x 10 3 . 

Figure [5] shows the results of the simulation, as a plot of system energy versus 
time. The mean energies of the system before and after the step change are 
indicated by a dashed horizontal line and an arrow, respectively, and the time 
of the change is indicated by the dashed vertical line. The simulation starts at 
time t = with the system at the initial mean energy estimate £\ = 10 4 . A 
total of n 1 = 3 x 10 4 flare events are simulated prior to the step change with 
the flare transition-rate coefficient 57oi, and Figure [5] shows the time history of 
the energy over this number of flares. The final jump transition occurs at time 
T = tx rj ni/Aoi = 3 x 10 6 . A total of 712 = 3 x 10 4 flares are also simulated after 
the step change during the interval of time £2 ~ ^2/^02 = 3 x 10 5 , and the time 
history of energy is again shown. The relaxation time estimate f \ = 1.1 X 10 3 is 
small compared with the scale of the time axis in Figure [5] Figure [5] illustrates 
the fluctuations of the system energy around the mean energy before and after 
the change, and the rapid relaxation of the system produced by a large increase 
in the rate. 

Figure |5] shows an expanded view in time of the evolution of the system energy 
during the interval of relaxation, for the simulation shown in Figure [5] The mean 
energies before and after the step change are shown by dashed horizontal lines, 
the time of the step change is shown by the dashed vertical line on the left, and 
the estimate f\ = 1.1 x 10 3 for the relaxation time is indicated by the interval 
between the dashed vertical lines. This figure shows the specific path E = E(t) 
to a new steady state taken in this simulation. 

Figure [7] presents the histograms of the frequency- energy distributions (upper 
row) and the waiting-time distributions (lower row) for the simulation shown in 
Figures [5] and [6] before (left column) and after (right column) the step change. 
The frequency-energy distributions arc plotted in a log-log representation, with 
the mean energy estimates £\ = 10 4 and £2 = 10 2 shown by the dashed vertical 
lines, and the simple power-law models for a system with a very large free energy 
[given by Equation (|37|) ] shown by the solid grey lines. The frequency-energy 
distributions follow the expected power laws for energies less than a rollover 
energy which is comparable to the mean energy estimate in each case. These 
results confirm the analytic arguments given in Section 12.31 The waiting-time 
histograms in the lower row of Figure [7J are shown together with the simple ex- 
ponential distributions corresponding to the estimates A01 = 0.01 and A02 = 0.1 
for the mean rates of events (solid grey lines). The results show that the waiting 
times follow approximate exponential (Poisson) distributions, as deduced by the 
analytic arguments in Section 12.31 However, there is a significant discrepancy 
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Figure 5. A Monte Carlo simulation of Case 1 showing the system energy versus time, for 
the choice of parameters ooi = 5 X 10 -3 and 5o2 = 5 X 10 -2 . A total of 3 X 10 4 events are 
simulated before and after the step change. The dashed horizontal line and the arrow show the 
mean energies before and after the change, respectively, and the time of the change is indicated 
by the dashed vertical line. 

between the slope of the simulation histogram and the slope of the simple model 
(corresponding to the rate estimate A02 = 0.1) after the change. In the simu- 
lation, the mean rates defined by the number of events divided by the elapsed 
time are A s i = ri\[i\ = 0.0099 (before the change) and A S 2 = n 2 /t2 = 0.090 
(after). The rate after the change is significantly less than the estimate given 
by Equation p9j) . which indicates that the approximation E 3> E c is not being 
met. Before the change, the mean energy is approximately £\ = 10 i E c , and 
Equation (j3"9")l provides a good approximation. After the change E\ = 10 2 E C , 
and the approximation is poorer. The rate is reduced because the low system 
energy prevents larger flares from occurring. However, the functional form of the 
waiting-time distribution remains approximately exponential. The lower row in 
Figure [7] also shows the simple exponential models corresponding to the mean 
rates A s i and X S 2 (dashed lines), and there is good agreement in both cases with 
the simulation histograms. 

Figure [5] illustrates part of the data from the same simulation, and is intended 
for qualitative comparison with the observational data for active region AR 
11029, shown in Figure The format of the figure is the same as Figure [2J The 
figure shows the events from the simulation with energy larger than E\ = 10 
which occurred in the interval of time between t = T — 30ri and t = T + 10ri, 
where T = 3.022 x 10 6 is the time of the change and f\ = 1.1 x 10 3 is the 
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Figure 6. An expanded view in time of the evolution of the system energy during the re- 
laxation interval, for the simulation shown in Figure [5] The estimate tj = 1.1 X 10 3 for the 
relaxation time is indicated. 

relaxation-time estimate. These choices are intended to crudely mimic observa- 
tional event selection above a threshold, and to include numbers of events before 
and after the change which provide good statistics but which are comparable 
with observational numbers for the soft X-ray events. The choices give Ni = 341 
events in total, with nn = 98 events before the change and nn = 243 after. The 
upper panel in Figure [5] shows the flare energies versus time, in the format of Fig- 
ure [H The middle panel shows the the flare cumulative number distribution for 
the selected events (grey diamonds) , which corresponds to the frequency- energy 
distribution according to Equation (U). The black line is the simple power-law 
model for a system with a very large free energy, corresponding to Equation (|37[) . 
Specifically, the black line is Cq(E) = N\{E / Ei)~~< +1 . The simulation events 
show a departure from the model Co(-Ei) at large energy due to the influence of 
the rollovers E\ and Ei (see the upper row of Figure [7]). The observed rollover for 
the selected events is qualitatively similar to the observations for active region 
AR 11029 (middle panel of Figure [5]). The lower panel in Figure |S] shows the 
cumulative waiting-time distribution for the selected events (grey diamonds), 
and the double exponential Poisson model (black curve) given by Equation Q : 



P(At) = -^-A lie xp(-AiiAi) + !^A 12 exp(-A 12 Ai), 



(44) 




where 



An 



30ti 



and A12 = 



10ti 



(45) 
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Figure 7. Histograms representing the frequency-energy (upper row) , and waiting-time (lower 
row) distributions before (left column) and after (right column) the step change, for the 
simulation shown in Figures [5] and [B] The dashed vertical lines in the upper row are the 
mean energy estimates and the solid grey curves are the power-law model for an active region 
with a a very large mean energy. The solid grey lines in the lower row are the the simple 
exponential (Poisson) models corresponding to the rates Aoi and A02, and the dashed lines are 
the exponential models with rates A s i and A s 2 (see text). 

are the rates of selected events before and after the change. This panel, which 
may be compared qualitatively with the lower panel in Figure [21 shows that the 
selected events follow the double power-law model. 

Figure[3]illustrates the relaxation process using an ensemble of 100 simulations 
each of tt.2 = 3 x 10 4 events generated for the time after the step change. 
In each simulation the system is started at time t = with energy equal to 
the estimate E\ = 10 for the mean energy of the system before the change, 
and flare events are simulated for the interval of time < t < 6ti. These 
choices simulate the relaxation of the system from the steady state existing 
before the change in the rate. The average energy over the ensemble of 100 
simulations is calculated for 600 equally spaced times during the simulation 
interval. Figure [9] shows the ensemble-mean energies versus time. The estimate 
Ti = 1.1 x 10 3 for the relaxation time is indicated by the dashed vertical line. The 
average energy decreases towards the mean energy estimate £2 = 10 2 (dashed 
horizontal line) following the change, reaching a statistically steady state after 
a time r s i m ~ 4000 — 5000 « 4ri. These results suggest that the estimate f\ 
underestimates the actual relaxation time, which is not surprising given the 
large increase in the flaring rate, and hence decrease in the mean energy of the 
system [see comments following Equation ([34]) ]. This figure also indicates that 
the estimate for the final mean energy is an underestimate. This is also not 
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Figure 8. Events in the simulation shown in Figurcs l5l7l largcr than size E\ = 10 which occur 
in the interval between 30 relaxation times before, and ten relaxation times after, the change 
in the rate. The format of the figure is the same as for the observational data for active region 
AR 11029 shown in Figure[2] 



surprising as it is an approximation, and in particular involves the assumption 
£2 > E c [see IWheatlandl (|2008p for the derivation] . 



4. Conclusions 

A stochastic description of the free energy of an active region, a 'jump-transition 
model' ( jWheatland and Glukhov, 1998| |Wheatland, 20081 |Wheatland, 2009] ) is 
applied to model an active region which exhibits time variation in its mean 
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Figure 9. An average of the energy of the system versus time over an ensemble of 100 
simulations constructed for the time after the step change. The simulations start at the mean 
energy estimate E\ = 10 4 for the time before the change. The estimate for the mean energy 
£2 = 10 2 after the change is shown by the solid horizontal line, and the dashed vertical line is 
the estimate t\ = 1.1 X 10 3 for the relaxation time. 

flaring rate. Time variation in flare productivity is commonly observed on the 
Sun, for example in NOAA solar active region AR 11029, which on October 26 
exhibited an approximate 10-fold increase in flaring, based on soft X-ray events 
compiled by the US Space Weather Prediction Center ([Wheatla nd, 2010) ). The 
general jump-transition model has previously been investigated only for time- 
independent situations, in which case it reproduces observed solar flare statistics 
(in particular the power-law flare frequency energy distribution, and a Poisson 
waiting-time distribution) . The interest here is with whether the time-dependent 
generalization of the model also succeeds in this respect. 

Time variation is incorporated in the model for two simple cases: 1. a step 
change in a coefficient ao(t) describing the rates of flaring; and 2. a step change 
in a coefficient /3o(i) describing the mean rate of energy of the system. Case 
1 may be appropriate to describe solar active region AR 11029. Analytic ar- 
guments are presented that predict the qualitative behavior of the model in 
response to the changes. In both cases the system adjusts by shifting to a new 
statistically stationary steady state following the change, over a relaxation time, 
which is analytically estimated. Steady states of the system have mean energies 
which may also be analytically estimated, as shown in I Wheatland a nd Glukhov 
(|1998[) . Provided the mean energy is sufficiently large, the steady-state system 
exhibits a power-law flare frequency energy distribution (for flare energies less 
than a rollover value corresponding approximately to the mean energy), and 
exhibits simple Poisson waiting-time statistics, as shown in [Wheatland (I2008P 
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and Wheatland (2009}. The analytic arguments given here predict that the 
same is true in time-dependent situations, although the waiting-time distribution 
corresponds to a time-dependent Poisson process with a rate determined by the 
model coefficient ao(t). For Case 1, the waiting-time distribution is in general 
a double exponential, corresponding to two intervals of flaring with a constant 
rate. 

Monte Carlo simulations of the stochastic differential equation formulation 
of the jump transition model (Wheatland, 2009) are presented for Case 1, for 
parameter choices intended to qualitatively model the soft X-ray event data 
for active region AR 11029. The simulations confirm the analytic arguments 
for the behavior of the model. The observed soft X-ray flare frequency-peak 
flux distribution for this active region showed evidence for an upper rollover 



( Wheatland, 2010 ), which was interpreted in terms of the finite amount of energy 
available for flaring in a small active region being depleted by an interval of 
very rapid flare production. This interpretation is consistent with the the jump- 
transition model, and the rollover is qualitatively reproduced in the simulations 
presented here. A double exponential waiting-time distribution is obtained in 
the simulations which mimics that observed for solar active region AR 11029. 

A quantitative comparison of data such as that for active region AR 11029 
with the jump transition model requires estimation of flare energies, and active 
region energies, from the data. However, soft X-ray peak fluxes are indicative 
only of flare energy, and we are currently unable to reliably estimate the magnetic 
free energy of a solar active region [see e.g. IPe Rosa et q"Tl(|2009[) for an attempt 
to obtain energy estimates for one active region via magnetic field modeling]. 
Nevertheless, in future work we will consider more detailed comparison of the 
model with observational data. 

An interesting aspect of the jump-transition model in this context is the model 
prediction of a departure from Poisson waiting-time statistics, if the system 
energy becomes sufficiently small, because large events are prevented from oc- 
curring, which reduces the overall mean rate of events. This can lead to departure 
from an exponential waiting-time distribution, even in the time-independent case 
(Wheatland, 2008[ [Wheatland, 2009[ ). However, the simulations presented here 
reveal a new effect. If the system mean energy has intermediate values, the 
overall rate of events is reduced by comparison with that for a system with 
a large mean energy, but the waiting-time distribution remains approximately 
exponential. Hence the waiting-time distribution may remain exponential even 
when there is detailed departure from Poisson statistics. This effect was not 



noticed in previous modeling (Wheatland, 2008 Wheatland, 20091, and also 



warrants further investigation. The soft X-ray events in AR 11029 appeared 
to follow Poisson occurrence in time based on agreement between the waiting- 



time distribution and the Poisson model ( Wheatland, 2010 1. However, the US 
SWPC event data suffer from significant event detection and selection problems 
due to the large and time varying soft X-ray background [discussed, for example, 
in I Wh eatland (20 01]) ] . Comparison with more sensitive data may be required to 
reveal detailed departure from Poisson occurrence. 

The advent of new high resolution and high cadence short-wavelength ob- 
servations of the Sun by the Advanced Imaging Assembly (AIA) on the Solar 
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Dynamics Observatory (SDO) should enable improved investigations of solar 
flare event statistics, and more careful comparison with flare statistics models 
such as the one in this paper. The results may provide new insights into the 
mechanisms of energy storage and release underlying the flare phenomenon. 

Acknowledgements The authors thank Kinwah Wu for comments on a draft of the paper. 
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